Different approaches to dimensionality reduction of DDM variables

source('~/Dropbox/PoldrackLab/SRO_DDM_Analyses/code/workspace_scripts/ez_fa_data.R')

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Warning: Unknown columns: `X`, `sub_id`, `subj_id`

Warning: Unknown columns: `X`, `sub_id`, `subj_id`

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Warning: Unknown columns: `X`, `sub_id`, `subj_id`

fig_path = '/Users/zeynepenkavi/Dropbox/PoldrackLab/SRO_DDM_Analyses/output/figures/'
source('~/Dropbox/PoldrackLab/SRO_Retest_Analyses/code/figure_scripts/figure_res_wrapper.R')
source('~/Dropbox/PoldrackLab/SRO_Retest_Analyses/code/helper_functions/make_rel_df.R')

To reduce the number of variables in a data-driven way (instead of just selecting the variables that went in to the ontology paper) I’ll apply the cleaning methods from the ontology pipeline: - dropping variables with r>0.85. - Remove outliers (>2.5 SD away) - transformation of non-normal variables (should be particularly useful for a set of variables with many response times) - Residualize Age and Sex effects

res_clean_test_data_ez_522_condition = res_clean_test_data_ez_522 %>%
  select((measure_labels %>% 
            filter(raw_fit == "EZ" & overall_difference == "condition"))$dv)

PCA on EZ variables of test data

Number of variables analysis

In this dataset there are more variables than observations.

Given the large number of measures some can possible be represented as linear combinations as others.

This would result in 0 or negative values in the eigenvalues of the correlation matrix and would mean that the correlation matrix is not positive definite.

It is hard to detect such multiple dependencies. Given the number of possible combinations of variables from all possible variables it is also impossible to compute the exact largest combination of variables that has a positive definite correlation matrix.

So instead I used a sampling approach. Below is a plot of the proportion of positive definite matrices out of 1000 samples of the given numbers of variables drawn from all variables.

The plot shows that there is a steady decrease in the proportion of positive definite matrices the more variables are used. It sharply drops to 0 from chance at 160 variables.

max_num_vars_summary = read.csv('/Users/zeynepenkavi/Dropbox/PoldrackLab/SRO_DDM_Analyses/input/max_num_vars_summary.csv')

p = max_num_vars_summary %>%
  ggplot(aes(as.factor(cur_num_vars), prop_det_pos))+
  geom_point()+
  geom_line(group=1)+
  xlab("Number of variables drawn")+
  ylab("Proportion of positive definite correlation matrices")+
  theme(axis.text.x = element_text(angle = 90))

ggsave('max_num_vars.jpeg',plot=p, device = 'jpeg', path = fig_path, width =12, height = 4, units = "in")

p

This means that any effort to reduce dimensionality using an eigenvalue decomposition will encounter problems.

For numerical efficiency and stability many PCA applications use singular value decomposition instead (e.g. see documentation for the princomp and prcomp functions in the base package).

To be consistent with the dimensionality reduction methods used in the ontology paper I will use the psych package functions for PCA and EFA. These rely on an eigenvalue decomposition but have a built-in function to smooth over negative eigenvalues that prints a warning but continues to execute and yield a solution.

But to confirm that the resulting factors are not problematic I will run an SVD and recover the same parameters as well.

Eigen on cor matrix

Eigenvalue decomposition on correlation matrix.

ez_t1_pca_3 = principal(res_clean_test_data_ez_522_condition, nfactors=3, rotate="oblimin")

Scree plot of first 10 components

data.frame(ez_t1_pca_3$values) %>%
  rename(eig = ez_t1_pca_3.values) %>%
  arrange(-eig) %>%
  mutate(var_pct = eig/sum(eig)*100,
         pc = 1:n()) %>%
  filter(pc<11)%>%
  ggplot(aes(factor(pc), var_pct))+
  geom_bar(stat="identity")+
  ylab("Percentage of variance explained")+
  xlab("Principal component")

Difference in percentage of variance explained is below 1% after the 3 component.

data.frame(ez_t1_pca_3$values) %>%
  rename(eig = ez_t1_pca_3.values) %>%
  arrange(-eig) %>%
  mutate(var_pct = eig/sum(eig)*100,
         pc = 1:n(),
         var_pct_shift = lead(var_pct),
         var_pct_diff = var_pct - var_pct_shift)

Visualizing the first three coponents and coloring them by the parameter type.

ez_t1_pca_3_loadings = as.data.frame(ez_t1_pca_3$loadings[])

ez_t1_pca_3_loadings[abs(ez_t1_pca_3_loadings)<0.3]=NA

tmp = ez_t1_pca_3_loadings %>%
  mutate(dv = row.names(.)) %>%
  select(dv, TC1, TC2, TC3) %>%
  mutate(num_loading = 3-(is.na(TC1)+is.na(TC2)+is.na(TC3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-TC1, -TC2, -TC3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "Drift rate",
                  ifelse(grepl("thresh", dv), "Threshold",
                         ifelse(grepl("non_dec", dv), "Non-decision", NA))),
         var_type = factor(var_type, levels=c("Drift rate", "Threshold", "Non-decision")))

p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank(),
        axis.text.y = element_blank())

ggsave('EZ_PCA_T1_3.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

p

SVD

Double checking the above factor solution using a singular value decomposition and oblimin rotation.

ez_svd <- svd(res_clean_test_data_ez_522_condition)
ncomp = 3
df <- nrow(res_clean_test_data_ez_522_condition) - 1
ez_t1_svd_rawLoadings = ez_svd$v[,1:ncomp] %*% diag(ez_svd$d/sqrt(df), ncomp, ncomp)
ez_t1_svd_rotatedLoadings <- GPArotation::oblimin(ez_t1_svd_rawLoadings)$loadings

ez_t1_svd_rotatedLoadings[abs(ez_t1_svd_rotatedLoadings)<0.3]=NA

tmp = data.frame(ez_t1_svd_rotatedLoadings) %>%
  mutate(dv = names(res_clean_test_data_ez_522_condition),
         X1 = -1*X1,
         X3= -1*X3) %>%
  select(dv, X1, X2, X3)%>%
  mutate(num_loading = 3-(is.na(X1)+is.na(X2)+is.na(X3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-X1, -X2, -X3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "Drift rate",
                           ifelse(grepl("thresh", dv), "Threshold",
                                  ifelse(grepl("non_dec", dv), "Non-decision", NA))),
                           var_type = factor(var_type, levels=c("Drift rate", "Threshold", "Non-decision")))

p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank(),
        axis.text.y = element_blank())

ggsave('EZ_SVD_T1_3.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

p

EFA on EZ variables of test data

How many factors should be extracted from the EFA? To answer this we run models extracting 2 to 50 components and rank them by BIC.

efa_ez_t1_comp_metrics = find_optimal_components(res_clean_test_data_ez_522_condition, fm = "minres", minc=2)

This suggests that a 13 component solution would be the best fitting model.

efa_ez_t1_comp_metrics

Fit the model suggested by the BIC comparison.

ez_t1_fa_13 = fa(res_clean_test_data_ez_522_condition, efa_ez_t1_comp_metrics$comp[1], rotate='oblimin', fm='minres', scores='tenBerge')

Fit the 3 factor model that is of theoretical interest.

ez_t1_fa_3 = fa(res_clean_test_data_ez_522_condition, 3, rotate='oblimin', fm='minres', scores='tenBerge')

Differences in BIC’s are difficult to interpret.

I tried to run a more formal model comparison between the 3 factor solution that we have theoretical reasons for and the 16 factor model that is selected based on BIC. Below is the result. I’m not quiet sure how to interpreted but it seems to suggest that the more complicated model is not significantly better.

anova(ez_t1_fa_3, ez_t1_fa_13)

So what would the 3 factor solution look like? Not surprisingly, very similar to the PC’s.

ez_t1_fa_3_loadings = as.data.frame(ez_t1_fa_3$loadings[])

ez_t1_fa_3_loadings[abs(ez_t1_fa_3_loadings)<0.3]=NA

tmp = ez_t1_fa_3_loadings %>%
  mutate(dv = row.names(.)) %>%
  select(dv, MR1, MR2, MR3) %>%
  mutate(num_loading = 3-(is.na(MR1)+is.na(MR2)+is.na(MR3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-MR1, -MR2, -MR3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "Drift rate",
                  ifelse(grepl("thresh", dv), "Threshold",
                         ifelse(grepl("non_dec", dv), "Non-decision", NA))),
                           var_type = factor(var_type, levels=c("Drift rate", "Threshold", "Non-decision")))


p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank())
        # axis.text.y = element_blank())

ggsave('EZ_FA_T1_3_w_labels.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 15, height = 15, units = "in")

p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank(),
        axis.text.y = element_blank())

ggsave('EZ_FA_T1_3.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

p

PCs vs Factors

Variance acconted for by PCA and EFA

ez_t1_pca_3$Vaccounted

                          TC1     TC2     TC3
SS loadings           15.9033 11.9575 9.58886
Proportion Var         0.1359  0.1022 0.08196
Cumulative Var         0.1359  0.2381 0.32008
Proportion Explained   0.4247  0.3193 0.25605
Cumulative Proportion  0.4247  0.7440 1.00000

ez_t1_fa_3$Vaccounted

                          MR1      MR2     MR3
SS loadings           15.3460 11.29760 8.86720
Proportion Var         0.1312  0.09656 0.07579
Cumulative Var         0.1312  0.22772 0.30351
Proportion Explained   0.4322  0.31815 0.24970
Cumulative Proportion  0.4322  0.75030 1.00000

Comparison of loadings

ez_t1_pca_3_loadings %>%
  mutate(dv = row.names(.)) %>%
  left_join(ez_t1_fa_3_loadings %>%
              mutate(dv = row.names(.)), by="dv") %>%
  gather(key, value, -dv) %>%
  mutate(factor_num = ifelse(grepl("1", key), 1, ifelse(grepl("2", key), 2, ifelse(grepl("3", key), 3, NA))),
         key = ifelse(grepl("TC", key), "PCA", ifelse(grepl("MR", key), "EFA", NA))) %>% spread(key, value) %>%
  ggplot(aes(PCA, EFA, col=factor(factor_num)))+
  geom_point()+
  geom_abline(aes(intercept=0, slope = 1), linetype=2)+
  facet_wrap(~factor_num)+
  theme(legend.position = "none",
        panel.grid = element_blank())

Artifical PCA

art_pca_scores = res_clean_test_data_ez_522_condition %>%
  mutate(sub_id = test_data_522$sub_id) %>%
  gather(dv, value, -sub_id) %>%
  left_join(measure_labels %>% select(dv, rt_acc), by="dv") %>%
  group_by(sub_id, rt_acc) %>%
  summarise(mean_val = mean(value)) %>%
  mutate(rt_acc = ifelse(rt_acc == "drift rate", "MR1", ifelse(rt_acc == "non-decision", "MR2", "MR3")))

Warning: attributes are not identical across measure variables;
they will be dropped

art_pca_scores

data.frame(ez_t1_fa_3$scores) %>%
   mutate(sub_id = test_data_522$sub_id) %>%
  gather(rt_acc, factor_score, -sub_id)%>%
  left_join(art_pca_scores, by=c("sub_id", "rt_acc")) %>%
  ggplot(aes(mean_val,factor_score, col=rt_acc))+
  geom_point()+
  geom_abline(aes(intercept=0, slope=1), linetype=2)+
  facet_wrap(~rt_acc, scales="free")+
  theme(legend.position = "none",
        panel.grid = element_blank())+
  xlab("Mean value of individual measures")+
  ylab("EFA")

Based on these results I ran 3 factor EFA’s on the remaining data.

EFA on HDDM variables of test data

Extract EZ variables

test_data_hddm = test_data %>%
  select(grep('hddm', names(test_data), value=T))

Remove variables that are correlated >0.85

clean_test_data_hddm = remove_correlated_task_variables(test_data_hddm)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Remove outliers (>2.5 SD away)

clean_test_data_hddm = as.data.frame(apply(clean_test_data_hddm, 2, remove_outliers))

Transform skewed variables

#Nothing to transform
numeric_cols = get_numeric_cols()
numeric_cols = numeric_cols[numeric_cols %in% names(clean_test_data_hddm) == T]
clean_test_data_hddm = transform_remove_skew(clean_test_data_hddm, numeric_cols)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 data positively skewed data were transformed:
dot_pattern_expectancy.hddm_thresh
stim_selective_stop_signal.hddm_thresh
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
0 positively skewed data could not be transformed successfully:

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Drop subject identifier column, mean impute and drop cols with no variance

clean_test_data_hddm_std = clean_test_data_hddm %>% mutate_if(is.numeric, scale)

#mean imputation
clean_test_data_hddm_std[is.na(clean_test_data_hddm_std)]=0

#drop cols with no variance
clean_test_data_hddm_std = clean_test_data_hddm_std %>%
  select_if(function(col) sd(col) != 0)

Residualize Age and Sex effects

clean_test_data_hddm_std = cbind(clean_test_data_hddm_std, demographics[,c("Age", "Sex")])

res_clean_test_data_hddm = residualize_baseline(clean_test_data_hddm_std)

Fit the 3 factor model that is of theoretical interest.

hddm_t1_fa_3 = fa(res_clean_test_data_hddm, 3, rotate='oblimin', fm='minres', scores='tenBerge')

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
= np.obs, : The estimated weights for the factor scores are probably
incorrect. Try a different factor extraction method.

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

The three factor model does not work as well for HDDM’s because of the 106 variables we have 75 are drift rates, 17 are threshold and 14 are non-decision times. This is because while all the parameters are fit separately for EZ only the drift rate parameter is allowed to vary by condition for the HDDM. So the model is mostly separating out drift rates from each other.

hddm_t1_fa_3_loadings = as.data.frame(hddm_t1_fa_3$loadings[])

hddm_t1_fa_3_loadings[abs(hddm_t1_fa_3_loadings)<0.3]=NA

tmp = hddm_t1_fa_3_loadings %>%
  mutate(dv = row.names(.)) %>%
  select(dv, MR1, MR2, MR3) %>%
  mutate(num_loading = 3-(is.na(MR1)+is.na(MR2)+is.na(MR3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-MR1, -MR2, -MR3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "drif rate",
                  ifelse(grepl("thresh", dv), "threshold",
                         ifelse(grepl("non_dec", dv), "non-decision", NA))))


p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank(),
        axis.text.y = element_blank())

ggsave('HDDM_T1_FA_3.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

p

Does it fit better if I use just the main drift rates and not the condition ones?

hddm_var_subset = c("adaptive_n_back.hddm_drift" , "attention_network_task.hddm_drift", "choice_reaction_time.hddm_drift", "directed_forgetting.hddm_drift", "dot_pattern_expectancy.hddm_drift" , "local_global_letter.hddm_drift", "motor_selective_stop_signal.hddm_drift",  "recent_probes.hddm_drift", "shape_matching.hddm_drift" , "simon.hddm_drift", "stim_selective_stop_signal.hddm_drift", "stop_signal.hddm_drift", "stroop.hddm_drift" , "threebytwo.hddm_drift", "adaptive_n_back.hddm_thresh", "attention_network_task.hddm_thresh", "choice_reaction_time.hddm_thresh", "directed_forgetting.hddm_thresh", "dot_pattern_expectancy.hddm_thresh.logTr", "local_global_letter.hddm_thresh", "motor_selective_stop_signal.hddm_thresh", "recent_probes.hddm_thresh", "shape_matching.hddm_thresh", "simon.hddm_thresh", "stim_selective_stop_signal.hddm_thresh.logTr", "stop_signal.hddm_thresh", "stop_signal.hddm_thresh_high", "stop_signal.hddm_thresh_low", "stop_signal.proactive_slowing_hddm_thresh", "stroop.hddm_thresh", "threebytwo.hddm_thresh", grep("non_dec", names(test_data_hddm), value = T))

res_clean_test_data_hddm_subset = res_clean_test_data_hddm %>% select(hddm_var_subset)

hddm_t1_fa_3_subset = fa(res_clean_test_data_hddm_subset, 3, rotate='oblimin', fm='minres', scores='tenBerge')

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

Then it works better but still not doing a great job separating out thresholds from non-decision times.

hddm_t1_fa_3_subset_loadings = as.data.frame(hddm_t1_fa_3_subset$loadings[])

hddm_t1_fa_3_subset_loadings[abs(hddm_t1_fa_3_subset_loadings)<0.3]=NA

tmp = hddm_t1_fa_3_subset_loadings %>%
  mutate(dv = row.names(.)) %>%
  select(dv, MR1, MR2, MR3) %>%
  mutate(num_loading = 3-(is.na(MR1)+is.na(MR2)+is.na(MR3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-MR1, -MR2, -MR3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "drif rate",
                  ifelse(grepl("thresh", dv), "threshold",
                         ifelse(grepl("non_dec", dv), "non-decision", NA))))


p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  ylab("Absolute Loading")+
  scale_color_identity()+
  theme(legend.position = "bottom",
        legend.title=element_blank())
        #axis.text.y = element_blank())

ggsave('HDDM_T1_FA_3_subset.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

knitr::include_graphics(paste0(fig_path, 'HDDM_T1_FA_3_subset.jpeg'))

EFA on EZ variables of retest data

Preparing the data

retest_data_ez = retest_data %>%
  select(grep('EZ', names(retest_data), value=T))

Remove variables that are correlated >0.85

clean_retest_data_ez = remove_correlated_task_variables(retest_data_ez)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Remove outliers (>2.5 SD away)

clean_retest_data_ez = as.data.frame(apply(clean_retest_data_ez, 2, remove_outliers))

Transform skewed variables

#Nothing to transform
numeric_cols = get_numeric_cols()
numeric_cols = numeric_cols[numeric_cols %in% names(clean_retest_data_ez) == T]
clean_retest_data_ez = transform_remove_skew(clean_retest_data_ez, numeric_cols)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Drop subject identifier column, mean impute and drop cols with no variance

clean_retest_data_ez_std = clean_retest_data_ez %>% mutate_if(is.numeric, scale)

#mean imputation
clean_retest_data_ez_std[is.na(clean_retest_data_ez_std)]=0

#drop cols with no variance
clean_retest_data_ez_std = clean_retest_data_ez_std %>%
  select_if(function(col) sd(col) != 0)

Residualize Age and Sex effects

clean_retest_data_ez_std = cbind(clean_retest_data_ez_std, demographics[,c("Age", "Sex")])

res_clean_retest_data_ez = residualize_baseline(clean_retest_data_ez_std)

There are two ways we can go about this:

1. Predicting factor scores for T2 data using the T1 factor solution
2. Fitting a new 3-factor model to the T2 data

Predict retest data using the 3 factor EFA from test data

Define helper functions to predict t2 factor scores using different scoring methods and calculating their ICCs with t1 factor scores.

predict_t2_fa_scores = function(t1_df = res_clean_test_data_ez, t2_df = res_clean_retest_data_ez, scores="tenBerge", nfactors = 3,rotate='oblimin',fm='minres', sub_ids = retest_data$sub_id){
  require(psych)
  require(tidyverse)
  t1_fa = fa(t1_df, nfactors, rotate=rotate, fm=fm, scores=scores)
  t1_fa_scores = data.frame(t1_fa$scores) %>%
    mutate(sub_id = sub_ids)

  t2_pred = predict(t1_fa, t2_df)
  t2_pred_scores = as.data.frame(t2_pred) %>%
    mutate(sub_id = sub_ids)
  return(list(t1_fa_scores=t1_fa_scores, t2_pred_scores=t2_pred_scores))
}

get_icc_for_score_type = function(scores){
  p = predict_t2_fa_scores(scores=scores)
  r_df = make_rel_df(p$t1_fa_scores, p$t2_pred_scores, metrics = c("icc2.1"))
  r_df = r_df %>%
    spread(dv,icc2.1) %>%
    mutate(scores = scores)
  return(r_df)
}

Using tenBerge our default option to calculate factor scores

ez_t2_fa_3_pred = predict_t2_fa_scores()

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

The determinant of the smoothed correlation was zero.
This means the objective function is not defined.
Chi square is based upon observed residuals.

The determinant of the smoothed correlation was zero.
This means the objective function is not defined for the null model either.
The Chi square is thus based upon observed correlations.

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
= np.obs, : The estimated weights for the factor scores are probably
incorrect. Try a different factor extraction method.

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs =
np.obs, : Matrix was not positive definite, smoothing was done

Plotting predicted tenBerge scores for T2 data using T1 FA against T1 FA scores.

ez_t2_fa_3_pred$t1_fa_scores %>%
  mutate(time = "t1_fa_scores") %>%
  rbind(ez_t2_fa_3_pred$t2_pred_scores %>% mutate(time="t2_pred_scores")) %>%
  gather(factor, score, -sub_id, -time) %>%
  spread(time, score) %>%
  ggplot(aes(t1_fa_scores, t2_pred_scores, col=factor))+
  geom_point()+
  geom_abline(aes(slope=1, intercept=0))+
  facet_wrap(~factor, scales='free')+
  theme(legend.position = "none")

No relationship! As a result ICC’s are 0 too.

BUT tenBerge is not the only way of calculating factor scores. psych offers different methods for calculating factor scores based on Grice (2001). Here I calculate the ICC’s of the three factors using other methods of calculating the scores.

Regardless of which method I use I get the following warning

The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.

I tried all five and only three gave me factor scores despite this warning. But they imply different reliabilities. Anderson and Harman yield high ICC’s for all three factors where tenBerge says it should be 0.

score_type_iccs = rbind(get_icc_for_score_type("tenBerge"),
      get_icc_for_score_type("Anderson"),
      get_icc_for_score_type("Harman"))

score_type_iccs

score_type_iccs %>%
  gather(factor, icc2.1, -scores) %>%
  ggplot(aes(scores, icc2.1, col=factor))+
  geom_point(size=2.5)+
  xlab('')+
  theme(legend.title=element_blank())

Plotting the scores using Anderson as an alternative example. Reflecting the high ICC’s now all factor scores for both time points are highly correlated.

ez_t2_fa_3_pred_bartlett = predict_t2_fa_scores(scores="Anderson")

ez_t2_fa_3_pred_bartlett$t1_fa_scores %>%
  mutate(time = "t1_fa_scores") %>%
  rbind(ez_t2_fa_3_pred_bartlett$t2_pred_scores %>% mutate(time="t2_pred_scores")) %>%
  gather(factor, score, -sub_id, -time) %>%
  spread(time, score) %>%
  ggplot(aes(t1_fa_scores, t2_pred_scores, col=factor))+
  geom_point()+
  geom_abline(aes(slope=1, intercept=0))+
  facet_wrap(~factor, scales='free')+
  theme(legend.position = "none")

Based on the warning in fitting the models I dug into the psych package to find the different weight matrix calculations depending on the score methods:

where x is data and r <- cor(x, use = “pairwise”) f <- loadings(factor_model) Phi <- factor_model$Phi

Thurston

S <- f %*% Phi
w <- try(solve(r, S), silent = TRUE)

tenBerge

L <- f %*% matSqrt(Phi)
r.5 <- invMatSqrt(r)
r <- cor.smooth(r)
inv.r <- try(solve(r), silent = TRUE)
C <- r.5 %*% L %*% invMatSqrt(t(L) %*% inv.r %*% L)
w <- r.5 %*% C %*% matSqrt(Phi)

Anderson

I <- diag(1, nf, nf)
h2 <- diag(f %*% Phi %*% t(f))
U2 <- 1 - h2
inv.U2 <- diag(1/U2)
w <- inv.U2 %*% f %*% invMatSqrt(t(f) %*% inv.U2 %*% r %*% inv.U2 %*% f)

Bartlett

I <- diag(1, nf, nf)
h2 <- diag(f %*% Phi %*% t(f))
U2 <- 1 - h2
inv.U2 <- diag(1/U2)
w <- inv.U2 %*% f %*% (solve(t(f) %*% inv.U2 %*% f))

Harman

m <- f %*% t(S)
diag(m) <- 1
inv.m <- solve(m)
w <- inv.m %*% f

I’m stuck here though and can’t figure out what the problem with the tenBerge weight matrix is or what to make of different reliability results depending on the scoring method.

Fit a 3 factor model separately for T2 data.

Instead of predicting T2 data from the T1 factor solution we can also fit a 3 factor model to the T2 data independently.

Conceptually this might be similar to fitting the DDM’s separately for two time points when generating parameter estimates instead of using the estimates from the first time point for each subject.

ez_t2_fa_3 = fa(res_clean_retest_data_ez, 3, rotate='oblimin', fm='minres', scores='tenBerge')

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(R): Matrix was not positive definite, smoothing was
done

Warning in cor.smooth(r): Matrix was not positive definite, smoothing was
done

The determinant of the smoothed correlation was zero.
This means the objective function is not defined.
Chi square is based upon observed residuals.

The determinant of the smoothed correlation was zero.
This means the objective function is not defined for the null model either.
The Chi square is thus based upon observed correlations.

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
= np.obs, : The estimated weights for the factor scores are probably
incorrect. Try a different factor extraction method.

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs =
np.obs, : Matrix was not positive definite, smoothing was done

Plot of independent 3 factor solution for T2 data looks very similar to that of T1 solution (with two factors switched names)

ez_t2_fa_3_loadings = as.data.frame(ez_t2_fa_3$loadings[])

ez_t2_fa_3_loadings[abs(ez_t2_fa_3_loadings)<0.3]=NA

tmp = ez_t2_fa_3_loadings %>%
  mutate(dv = row.names(.)) %>%
  select(dv, MR1, MR2, MR3) %>%
  mutate(num_loading = 3-(is.na(MR1)+is.na(MR2)+is.na(MR3) ) ) %>%
  filter(num_loading!=0) %>%
  select(-num_loading) %>%
  arrange(-MR1, -MR2, -MR3) %>%
  mutate(order_num = 1:n(),
         dv = reorder(dv, -order_num)) %>%
  select(-order_num) %>%
  gather(Factor, Loading, -dv) %>%
  na.exclude() %>%
  mutate(neg_load = factor(ifelse(Loading>0,"NA","#000000")),
         var_type = ifelse(grepl("drift", dv), "drif rate",
                  ifelse(grepl("thresh", dv), "threshold",
                         ifelse(grepl("non_dec", dv), "non-decision", NA))))


p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  scale_color_identity()+
  ylab("Absolute Loading")+
  theme(legend.position = "bottom",
        legend.title=element_blank())

ggsave('EZ_FA_T2_3_w_labels.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 15, height = 15, units = "in")

p = tmp%>%
  ggplot(aes(dv, abs(Loading), fill=var_type, col = neg_load))+
  geom_bar(stat = "identity", alpha=0.5)+
  facet_wrap(~Factor, nrow=1)+
  coord_flip()+
  xlab("")+
  scale_color_identity()+
  ylab("Absolute Loading")+
  theme(legend.position = "bottom",
        legend.title=element_blank(),
        axis.text.y = element_blank())

ggsave('EZ_FA_T2_3.jpeg',plot=p, device = 'jpeg', path = fig_path, width = 10, height = 12, units = "in")

p

T1_fit scores vs T2_fit scores

ez_t2_fa_3_scores = data.frame(ez_t2_fa_3$scores)
ez_t2_fa_3_scores$sub_id = retest_data$sub_id

data.frame(ez_t1_fa_3$scores) %>%
  mutate(sub_id = test_data_522$sub_id) %>%
  left_join(ez_t2_fa_3_scores, by = "sub_id") %>%
  gather(key, value, -sub_id) %>%
  separate(key, c("par", "model"), sep = "\\.") %>%
  mutate(model = ifelse(model=="y", "T2_fit","T1_fit")) %>%
  spread(model, value)%>%
  ggplot(aes(T1_fit, T2_fit))+
  geom_point()+
  facet_wrap(~par, scales= "free")+
  geom_abline(aes(slope=1, intercept=0))

Warning: Removed 1116 rows containing missing values (geom_point).

T1_fit vs T2_fit loadings

ez_t1_fa_3_loadings = data.frame(ez_t1_fa_3$loadings[]) %>%
  rename(drift_rates = MR1, thresholds = MR2, non_decs = MR3) %>%
  mutate(dv = row.names(.))

ez_t2_fa_3_loadings = data.frame(ez_t2_fa_3$loadings[]) %>%
  rename(drift_rates = MR1, thresholds = MR3, non_decs = MR2)%>%
  mutate(dv = row.names(.))

ez_t1_fa_3_loadings %>%
  left_join(ez_t2_fa_3_loadings, by = "dv") %>%
  gather(key, value, -dv) %>%
  separate(key, c("par", "model"), sep = "\\.") %>%
  mutate(model = ifelse(model == "x", "T1_fit", "T2_fit"))%>%
  spread(model, value) %>%
  ggplot(aes(T1_fit, T2_fit))+
  geom_point()+
  geom_abline(aes(slope=1, intercept=0))+
  facet_wrap(~par)

Model time as a latent variable

Dave’s response to my confusions on this:

The standard way to evaluate test-retest correlations with latent variables is to fit a latent variable for each time and then allow the latent variables to correlate. That is the test-retest correlation for the best, “true” score measure at each time point.

Based on this I tried modeling all the data together using timepoint it was collected in as a latent variable. The details of those analyses are here

EFA on HDDM variables of retest data

retest_data_hddm = retest_data %>%
  select(grep('hddm', names(retest_data), value=T))

Remove variables that are correlated >0.85

clean_retest_data_hddm = remove_correlated_task_variables(retest_data_hddm)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Remove outliers (>2.5 SD away)

clean_retest_data_hddm = as.data.frame(apply(clean_retest_data_hddm, 2, remove_outliers))

Transform skewed variables as in T1

clean_retest_data_hddm = clean_retest_data_hddm %>%
  mutate(dot_pattern_expectancy.hddm_thresh.logTr = log(dot_pattern_expectancy.hddm_thresh),
stim_selective_stop_signal.hddm_thresh.logTr = log(stim_selective_stop_signal.hddm_thresh)) %>%
  select(-dot_pattern_expectancy.hddm_thresh,-stim_selective_stop_signal.hddm_thresh)

Drop subject identifier column, mean impute and drop cols with no variance

clean_retest_data_hddm_std = clean_retest_data_hddm %>% mutate_if(is.numeric, scale)

#mean imputation
clean_retest_data_hddm_std[is.na(clean_retest_data_hddm_std)]=0

#drop cols with no variance
clean_retest_data_hddm_std = clean_retest_data_hddm_std %>%
  select_if(function(col) sd(col) != 0)

Residualize Age and Sex effects

clean_retest_data_hddm_std = cbind(clean_retest_data_hddm_std, demographics[,c("Age", "Sex")])

res_clean_retest_data_hddm = residualize_baseline(clean_retest_data_hddm_std)

Select only the vars that went in to the limited set of hddm variables

res_clean_retest_data_hddm_subset = res_clean_retest_data_hddm %>%
  select(hddm_var_subset)

Predict using the 3 factor EFA from T1

hddm_t2_fa_3_subset = predict(hddm_t1_fa_3_subset, res_clean_retest_data_hddm_subset)